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Abstract 

We present a first-principles method for calculating the charging energy of a molecular 
single-electron transistor operating in the Coulomb blockade regime. The properties of the 
molecule are modelled using density-functional theory, the environment is described by a con- 
tinuum model, and the interaction between the molecule and the environment are included 
through the Poisson equation. The model is used to calculate the charge stability diagrams of 
a benzene and Ceo molecular single-electron transistor. 



Introduction 

The use of non-equilibrium Greens functions (NEGF) in connection with density-functional theory 
(DFT)™ 

or semi-empirical modelsl^^ has been highly successful in modelling coherent trans- 
port in various types of molecular junctions. However, in the case of molecular single-electron 
transistors (SET), the transport is incoherent,!^ and another approach is needed. Kaasbjerg et al. 
introduced a semi-empirical model for simulating the properties of molecular SETs. In particular 
*To whom correspondence should be addressed 
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they showed the importance of including renormalization of the molecular charge states due to the 
polarization of the environment. 

In this paper we extend this framework to be included within a density-functional theory de- 
scription of molecular SET's operating in the coulumb blockade regime. We use the model to 
calculate the charging energy of benzene and €50 in an electrostatic environment resembling a 
molecular SET geometry. We calculate the charging energy as function of an external gate poten- 
tial, and from this we obtain the charge stability diagram of the two respective molecules. 

The outline of the paper is as follows. In the first section we describe the basics of an SET, and 
in the following section the DFT framework for modelling the device. We next present calculations 
of the properties of benzene and C(,o in an electrostatic environment, and in the final section we 
summarize the results. 



Basic theory of a molecular single electron transistor 



[figure] [1] [] 1 1 schematically illustrates the geometry of a nanoscale molecular transistor. The 
geometry consists of metallic source and drain electrodes, and a molecular island coupled with the 
two electrodes. Electrons can propagate from source to drain through the island. 

If the island is strongly coupled with the source and drain electrodes, the electrons will stay 
a very short time on the island, and cannot localize but will move coherently through the system. 
This is the regime where we can use the coherent transport model for simulating the electrical 



properties of the system. This is the situation shown in [figure] [1][]1 d, and we note that a current 
can flow through the system even when the island does not have any electronic states within the bias 
window, as illustrated in the figure by the finite lifetime-broadened lowest unoccupied molecular 
orbital (LUMO), which the carrier can use for propagating from left to right. 

In this paper, we will investigate the regime where the island is weakly coupled with the elec- 
trodes. In this case, the electron tunnels from the source to the island and stays there for sufficiently 
long time to localize. The electron thereby loses all information about its original quantum state. 
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b) Coherent Tunneling Regime 
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c) Sequential Tunneling Regime 
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Figure 1: (Color online) a) Schematics of a nanoscale transistor. Electrons propagate from the 
source to the drain through an island. The energies of the electronic states on the island can 
be controlled by an electrostatic gate, b) In the strong coupling regime, the electron propagates 
coherently through the device, and the electronic states of the island have a finite lifetime and 
therefore a finite width, c) In the weak coupling regime, the electron localizes on the island and 
propagates by sequential tunneling. The electronic states on the island are in this case discrete. 
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Thus, the subsequent tunneling process from the island to the drain electrode will be independent 
from the tunneling process into the island. This transport mechanism is referred to as sequential 



tunneling and is depicted in [figure] [1][]1 Electron transport is in this case only possible if the is- 
land has an electronic level within the bias window, as illustrated by the electron affinity level (EA) 
in 



[figure] [!][]![ As in the coherent case, the position of the molecular levels and in particular the 
electron affinity level can be adjusted by an external gate potential, and by appropriate tuning, the 
island can in this regime thus be opened or closed for transport. For source-drain voltages below 
the charging energy of the island, there will only be one energy level within the bias window, and 
the system will work as a single electron transistor, as desired for our present discussion. 

The energy balance in the weak coupling regime 

In the following we will focus on the weak coupling regime where the transport is described by 
sequential tunneling. We introduce the function E^^^^^'^{N), which gives the total energy of the 
island as function of the number of electrons on the island. We also introduce similar energy 
functions for the source and drain electrodes, E^°^'^^^{N) and E'^^'^^{N). 

For the electron to move from the source electrode onto the island, the electron must have a 
lower energy on the island, i.e. 

^source ^^islandj^^^) > £Source^^_ j ) ^ ^island^^ ^ 

where and M are the initial number of electrons on the island and in the source electrode, 
respectively. 

Moreover, in order to move the electron from the island to the drain electrode, it must have a 
lower energy in the drain electrode, 

^drain ^ ^island (A^ + 1 ) > ^drain + I) + E^^^^"'* (A^) , (2) 

where K is the initial number of electrons in the drain electrode. 
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The maximum energy of the electron in the source electrode is —W + eV /2, where W is the 
work function of the electrode and V the applied bias. Assuming that the electron with maximum 
energy tunnels onto the island, then we have 

^source _ ^source _ j ^ _ _^ + gy/2. (3) 

Using the above tunneling criterion, we obtain the condition 

-W + eV/2 + E''^'""^{N)>E''^'"''^{N+l). (4) 

Similarly, —W — eV /lis the minimum energy of an electron in the drain electrode, and thus 

^island 1) > -W-eV/2 + E''^'"'\N). (5) 

The requirement for a current to flow in the device is therefore 

e\V\/2 > A£^^>^"'^(A^) > -e\V\/2, (6) 

where AE''^'''"^{N) = E''^'^'^{N+ 1) -E'''^'^^{N) is the charging energy of the island. 

In the following we will calculate the charging energy of two different molecular SETs and 



use 



[equation] [6] [] 6 to obtain the so-called charge stability diagram, which shows the number of 



charge states inside the bias window as function of the gate and source-drain voltages. 
Total energy of a molecule in an electrostatic environment 

In this section we will discuss the calculation of the total energy of a molecule in an electrostatic 
environment. All calculations are performed using the commercial software package Atomistix 
ToolKit (ATK).'^The DFT model in ATK is based on pseudopotentials with numerical localized 
basis functions using the method outlined by Soler et al. In this framework, a compensation 
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charge pf°"^^{r) is introduced for each atomic site. The compensation charge has the same charge 
Zi as the pseudopotential, and is used to screen the electrostatic interactions. 
We now introduce the electron difference density 

dn{r)=nir)-l^pr'ir): (7) 

i 

where n{r) is the total charge density of the system. We also introduce the screened ("neutral 
atom") local pseudopotential 

^NA/„x _ T/loc/„\ / Pi ) j„/ 



V;^A(r) = Vl°'{r) - / -jv^dr', (8) 



r — r 



where V'/°'^ is the local pseudopotential at site /. 

Following Ref. [151 we rearrange the terms to obtain the DFT total-energy functional 



E[n] = r[n] + i 1 5y«(r)5n(r)Jr + J '£v!'\r)dn{r)dr+^-'£u,j, (9) 

where T is the one-electron kinetic energy, E^'^ the exchange-correlation energy, and the three last 
terms are the rearranged electrostatic terms. The first of these is the relative Hartree energy, and is 
obtained from the difference Hartree potential dV^{r) which is calculated by solving the Poisson 
equation for the difference density dn{r). The next term describes the relative external energy, 
while the last term contains all the electrostatic interactions which do not depend on the electron 
density. The last term is calculated from Vj°^, pf°^^, and Z,-, and since the Poisson equation is 
hnear, this term can be decomposed into the pair-potential Uij. 

We next extend the total-energy functional to include interactions with a number of dielectric 



and metallic regions surrounding the molecular system, [figure] [2] []2 illustrates a typical molecu- 
lar single electron transistor geometry where a benzene molecule is positioned on top of a dielectric 
material and surrounded by three metallic electrodes. Within the metallic regions the potential is 
fixed to the applied voltage on each respective electrode. 
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Figure 2: (Color online) A benzene molecule in the SET environment considered in this paper. The 
electrostatic environment models two metallic electrodes (source and drain) on top of a dielectric 
substrate with a metallic back-gate. The contour plot shows the induced electrostatic potential for 
a gate voltage of 2 V and zero source-drain bias. 
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Solving the Poisson equation without the molecule present, we obtain the external potential 
from the electrostatic environment, 

-V- [e(r)Vy'='''(r)] =0. (10) 

Adding the molecule to the geometry, we again solve the Poisson equation, but now with the 
electron difference density on the molecule, to obtain the total difference Hartree potential 

-V - [e(r) V5y"+'='^'(r)] = 5n{r). (11) 

Finally, we define the molecular part of the total difference Hartree potential 

5V%) = 5y (r) - V^'^'ir) , (12) 



and this is the difference Hartree potential which enters into [equation] [9][]9[ Following Neuge 



bauer and Scheffler,''^we add the energy contribution from the external field through the term 

AE= [ y-'(r)n(r)J(r) -£y-'(R,-)Z,-, (13) 
where R/ is the position of site /, and the valency of the pseudopotential at site /. 



Thus, the total energy is given by adding the contributions from [equation] [9] [] 9 and [equa 



tion][13][]13 Note that this is only exact if the last electro- static pair-potential term in [equa 



tion][9][]9 is unaffected by the electrostatic surroundings, which is the case if the compensation 
charge and the screened local pseudopotential do not overlap with the metallic and dielectric re- 
gions. 
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Results 



We will now present results for the charging energy of a benzene and Cso molecule in an SET 
geometry. To obtain some reference energies, we will first calculate the charging energies of the 
isolated molecules. We obtain the charging energy by performing self-consistent calculations for 
the and 1 charge states of the isolated molecule and subtracting their total energies. For 
the calculation we use non-polarized DFT in the local density approximation (LDA)'i^and expand 
the wavefunctions in a double-zeta polarized (DZP) basis set. The molecular geometries where 
obtained by relaxing the molecules in the neutral state. 

Table 1 : Experimental ionization (I) and affinity (A) energies for benzene and Cgo molecules com- 
pared with theoretical DFT-LDA values. The theoretical values denoted "isolated" are obtained 
from the total energies of the charged isolated molecule. The "SET" values are obtained by cal- 
culating the total energy of the charged molecule in the electrostatic surrounding illustrated in 



[figure] [2] []2 at zero gate and source-drain bias. 



benzene 




Ei 






Exp. 




9.25IMI 






isolated 


15.73 


9.15 


-2.34 


-8.39 


SET 


7.70 


5.41 


-2.26 


-4.88 


Ceo 




El 


Ea 


E-^ 

^A 


Exp. 


11.33^'^ 




2.65^ 




isolated 


10.09 


6.84 


1.90 


-1.32 


SET 


7.24 


5.89 


2.85 


1.53 



[table] [1][]1 1 shows the calculated charging energies. For benzene there is excellent agreement 



with the experimental results, while the results for Ceo are about 1 eV too small. To investigate 
possible origins of this discrepancy we performed calculations where Cgo was allowed to spin 
polarize and relaxed the molecule also in the charged state. We found that such effects change the 
total energy by less than 0. 1 eV, and can therefore safely be disregarded. 

We next set up the molecules in the SET electrostatic environment. The geometry is illustrated 

m 



[figure] [2] []2 It consists of a metallic back-gate, and above the gate there is 3.8 A of dielectric 
material with dielectric constant IOCq. The molecule is positioned 1.2 A above the dielectric. To 
the left and right of the molecule are metallic source-drain electrodes, and the distance between 
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the molecule and the electrodes is 2.8 A. We note that for a typical metal surface the image plane 
is 2 A above the surface,!^ and to compare with atomic adsorption geometries this length must be 
added to the above distances. 



[table] [1][]1 lists the charging energies of the molecules in the SET environment with zero 
potential at the gate electrode. We see that for benzene the charging energy is strongly reduced, 
while the effect is smaller for Cgo- The reduction in charging energy arises from the screening of 
the charged molecule by the surrounding dielectric and metal electrodes. I^^l 

For benzene we may compare with the GW calculations by Neaton et al.^ They find the gas 
phase value for ^lumo — ^^homo = 10.51 eV and for benzene adsorbed on graphene ^lumo — 
£^HOMO = 7.35 eV, which is comparable with our values of Ei — Ea = 1 1.49 eV in the gas phase 
and El — E^ = 7.67 eV in the SET environment. 

We next calculate the total energy of the different charge states of the SET system as function of 



the gate potential. The results are shown in [figure] [3][]3 The total energy includes the reservoir 



energy qW, where q is the charge of the molecule and W is the work function of the electrode; we 
use the value W = 5.28 eV which models a gold electrodcl^S 



[figure] [3] []3 shows that the neutral molecule has the lowest energy at zero gate potential. At 
negative gate potentials the positive charge states are stabilized, while the negative charge states 
are stabilized at positive gate potentials. This is in agreement with that HOMO and LUMO levels 
follow —eVc, thus, at positive bias the LUMO level gets below the electrode Fermi level and 
attracts an electron, and the molecule becomes negatively charged. At negative gate potentials the 
HOMO level gets above the electrode Fermi level and an electron is escaping from the molecule, 
which becomes positively charged. 

To understand the dependence between the total energy and the gate potential we fit a quadratic 
function to the data 

E = aqVG + l5{eVGf. (14) 

Note that we assume the linear term to be proportional to the charge q on the molecule, while 
the quadratic term arises from polarization of the molecule and therefore is independent of q. By 
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Figure 3: (Color online) a) Total energy as function of the gate voltage for the benzene molecule, 
in the SET environment. Different curves are for charge states —2, —1,0, 1,2. c) Charge stability 
diagram; the color shows the number of charge states within the bias window for a given gate 
voltage and source-drain bias, dark blue(O), blue(l), green(2), yellow(3), red(4). b) and d) shows 
similar plots for the C(,q molecule, in the SET environment. 
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fitting the data in [figure] [3][]3 we find for benzene a = 0.62 eV, /3 = —0.003 eV \ and for €50 
a = 0.38 eV, /3 = —0.025 eV ^ where the variation with the charge state is ~ 0.01 eV for a and 
-0.001 eV-1 for/3. 

Thus, benzene is stronger coupled with the gate than Cgo, because the benzene atoms on aver- 
age are closer to the dielectric substrate. Therefore, the benzene molecule shows an almost linear 
relationship between the total energy and the gate potential, since all atoms are almost identically 
shifted by the gate potential. 

For Cgo, on the other hand, the relationship between the total energy and the gate potential is 
non-linear. For this molecule the atoms closest to the dielectric region screen the gate potential for 
the rest of the molecule, and thus the gate coupling becomes smaller. A difference in the charges 
on different atoms in the molecule gives rise to a molecular dipole, and it is this polarization energy 
that gives the second-order contribution to the total energy. 



From the total energies we can finally calculate the charge stability diagram using [equa- 



tion] [6] []6 The result is show in [figure][3][]3| :,d. The different colors show the number of charge 



states in the bias window. We see that the excitation energy for Ceo is smaller than for benzene. 
For both systems the excitation energy of the second electron is much smaller than for the first 
electron. We also note that the non-linear dependence of the total energy on the gate potential for 
C60 is not observable in the charge stability diagram. This is because the charge stability diagram 
only depends on the energy differences between the charge states, and the second-order term in 



[equation] [1 4] [] 14 is independent on the charge state. 



Conclusions 

We have in this paper demonstrated the use of density-functional theory for calculating the charg- 
ing energy of a molecule in a metallic environment that models the geometry of a molecular 
single-electron transistor. We find that the metallic environment reduces the charging energy of 
the molecules, in agreement with and Huckel-^ calculations. We calculated the charging 
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energy as function of the gate potential and from this obtained the charge stability diagram. The 
simulations show how DFT can be used to gain new insight into the properties of molecular single- 
electron transistors operating in the coulumb blockade regime. 
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